Exploring PACE Surface Reflectance Data

1. PACE OCI L2 Surface Reflectance

Summary

Learning Objectives

Contents

  • Read PACE L2 Surface Reflectance data
  • Mask clouds and water from images
  • Display spectral signatures for selected locations
  • Calculate unique spectral indices

1. Setup

# !conda install cf_xarray -y
# Import required libraries
import cartopy
import cf_xarray
import earthaccess
import numpy as np
import rasterio as rio
import hvplot.xarray  
import matplotlib.pyplot as plt
import xarray as xr
import pyproj
import holoviews as hv
from io import BytesIO
auth = earthaccess.login(persist=True)
tspan = ("2025-03-01", "2025-03-20")
bbox = (-86.512, 31.70, -86.512, 34.938)
 
results = earthaccess.search_data(
    short_name="PACE_OCI_L2_SFREFL",
    temporal=tspan,
    bounding_box=bbox,
)
links = [result.data_links()[0] for result in results]
fs = earthaccess.get_fsspec_https_session()

To open the PACE L2 Surface Reflectance data, we will use a combination of fsspec and BytesIO to read the data directly into memory. These files are roughly 750 mb each. The fsspec https session authenticates via NASA Earthdata credentials, allowing access to the data and BytesIO allows us to read the data into memory without saving it to disk.


with fs.open(links[0]) as f:
    file_content = BytesIO(f.read())

These files are heirarchichal netCDF files, so we will use the datatree function to read in and view all of the groups. Using the datatree representation, we can browse through the groups and variables in the file, which is useful for understanding the structure of the data, and for selecting what we want to work with.

datatree = xr.open_datatree(file_content, decode_timedelta=False)
datatree
<xarray.DatasetView> Size: 0B
Dimensions:  ()
Data variables:
    *empty*
Attributes: (12/45)
    title:                             OCI Level-2 Data SFREFL
    product_name:                      PACE_OCI.20250301T183734.L2.SFREFL.V3_...
    processing_version:                3.0
    history:                           l2gen par=/data11/sdpsoper/vdc/vpu10/w...
    instrument:                        OCI
    platform:                          PACE
    ...                                ...
    geospatial_lon_max:                -74.50793
    geospatial_lon_min:                -105.216805
    startDirection:                    Ascending
    endDirection:                      Ascending
    day_night_flag:                    Day
    earth_sun_distance_correction:     1.0183274745941162

TODO Explanation of data structure and selection of these variables

ds = xr.merge(
    (
        datatree.ds,
        datatree["geophysical_data"].ds[["rhos", "l2_flags"]],
        datatree["sensor_band_parameters"].coords,
        datatree["navigation_data"].ds.set_coords(("longitude", "latitude")).coords,
    )
)
ds
<xarray.Dataset> Size: 1GB
Dimensions:        (number_of_lines: 1709, pixels_per_line: 1272,
                    wavelength_3d: 122)
Coordinates:
  * wavelength_3d  (wavelength_3d) float64 976B 346.0 351.0 ... 2.258e+03
    longitude      (number_of_lines, pixels_per_line) float32 9MB ...
    latitude       (number_of_lines, pixels_per_line) float32 9MB ...
Dimensions without coordinates: number_of_lines, pixels_per_line
Data variables:
    rhos           (number_of_lines, pixels_per_line, wavelength_3d) float32 1GB ...
    l2_flags       (number_of_lines, pixels_per_line) int32 9MB ...
Attributes: (12/45)
    title:                             OCI Level-2 Data SFREFL
    product_name:                      PACE_OCI.20250301T183734.L2.SFREFL.V3_...
    processing_version:                3.0
    history:                           l2gen par=/data11/sdpsoper/vdc/vpu10/w...
    instrument:                        OCI
    platform:                          PACE
    ...                                ...
    geospatial_lon_max:                -74.50793
    geospatial_lon_min:                -105.216805
    startDirection:                    Ascending
    endDirection:                      Ascending
    day_night_flag:                    Day
    earth_sun_distance_correction:     1.0183274745941162
# def tap_bounds(to_tap, res, direction):
#     """Tap bounds to regular grid."""
#     mult = 1 if to_tap >= 0 else -1
#     abs_to = abs(to_tap)
#     abs_res = abs(res)
#     mod_val = abs_to % abs_res
#     if mod_val == 0:
#         adj = 0
#     elif direction == "up" and mult == 1:
#         adj = abs_res - mod_val
#     elif direction == "down" and mult == 1:
#         adj = -mod_val
#     elif direction == "up" and mult == -1:
#         adj = -mod_val
#     elif direction == "down" and mult == -1:
#         adj = abs_res - mod_val
#     else:
#         raise ValueError("direction must be 'up' or 'down'")
#     return mult * (abs_to + adj)

# def build_pace_glt(ds, target_resolution=(0.01,-0.01), radius_of_influence=5000):
#     """
#     Build GLT array from xarray Dataset of PACE navigation data.

#     Parameters:
#     - ds: xarray.Dataset with 'latitude' and 'longitude' variables.
#     - target_resolution: tuple (res_x, res_y) in degrees; res_y can be negative.
#     - target_extent_ul_lr: optional tuple (min_x, max_y, max_x, min_y); else computed from data.
#     - radius_of_influence: float, max distance in meters for nearest neighbor.

#     Returns:
#     - np.ndarray of shape (2, y_size, x_size): [0,:,:] = column indices (1-based), [1,:,:] = row indices (1-based); -9999 for invalid.
#     """
#     lats = ds['latitude'].values
#     lons = ds['longitude'].values
#     lines, pixels = lons.shape

#     min_x = np.nanmin(lons)
#     max_y = np.nanmax(lats)
#     max_x = np.nanmax(lons)
#     min_y = np.nanmin(lats)

#     res_x, res_y = target_resolution

#     # Tap bounds for alignment
#     min_x = tap_bounds(min_x, res_x, "down")
#     max_y = tap_bounds(max_y, res_y, "up")
#     max_x = tap_bounds(max_x, res_x, "up")
#     min_y = tap_bounds(min_y, res_y, "down")

#     x_size = int(np.ceil((max_x - min_x) / abs(res_x)))
#     y_size = int(np.ceil((max_y - min_y) / abs(res_y)))

#     # Define swath and area
#     swath_def = pyresample.geometry.SwathDefinition(lons=lons, lats=lats)
#     area_extent = (min_x, min_y, max_x, max_y)
#     area_def = pyresample.geometry.AreaDefinition(
#         'glt_grid', 'GLT grid', 'latlong',
#         {'proj': 'latlong', 'ellps': 'WGS84'},
#         x_size, y_size, area_extent
#     )

#     # Get neighbour info using pyresample's KDTree (geospatial-aware)
#     valid_input, valid_output, index_array, distance_array = pyresample.kd_tree.get_neighbour_info(
#         swath_def, area_def, radius_of_influence=5000, neighbours=1
#     )

#     # Set invalid to -1
#     index_array = index_array.flatten().astype(np.int64)  # Ensure flat
#     index_array[~valid_output.flatten()] = -1

#     # Unravel to col (pixel), row (line), 1-based
#     glt_col = np.full(index_array.shape, -9999, dtype=int)
#     glt_row = np.full(index_array.shape, -9999, dtype=int)
#     valid_map = index_array != -1
#     glt_row[valid_map] = (index_array[valid_map] // pixels) + 1
#     glt_col[valid_map] = (index_array[valid_map] % pixels) + 1

#     # Reshape and stack
#     glt = np.stack((glt_col.reshape(y_size, x_size), glt_row.reshape(y_size, x_size)))

#     return glt
#Check the wavelengths available for PACE 
ds["wavelength_3d"]
<xarray.DataArray 'wavelength_3d' (wavelength_3d: 122)> Size: 976B
array([ 346.,  351.,  356.,  361.,  366.,  371.,  375.,  380.,  385.,  390.,
        395.,  400.,  405.,  410.,  415.,  420.,  425.,  430.,  435.,  440.,
        445.,  450.,  455.,  460.,  465.,  470.,  475.,  480.,  485.,  490.,
        495.,  500.,  505.,  510.,  515.,  520.,  525.,  530.,  535.,  540.,
        545.,  550.,  555.,  560.,  565.,  570.,  575.,  580.,  586.,  615.,
        620.,  625.,  630.,  635.,  640.,  642.,  645.,  647.,  650.,  652.,
        655.,  657.,  660.,  662.,  665.,  667.,  670.,  672.,  675.,  677.,
        679.,  682.,  697.,  699.,  702.,  704.,  707.,  709.,  712.,  714.,
        719.,  724.,  729.,  734.,  739.,  742.,  744.,  747.,  749.,  752.,
        754.,  772.,  774.,  779.,  784.,  789.,  794.,  799.,  804.,  809.,
        814.,  819.,  824.,  829.,  835.,  840.,  845.,  850.,  855.,  860.,
        865.,  870.,  875.,  880.,  885.,  890.,  895., 1038., 1249., 1618.,
       2131., 2258.])
Coordinates:
  * wavelength_3d  (wavelength_3d) float64 976B 346.0 351.0 ... 2.258e+03
Attributes:
    long_name:  wavelengths
    units:      nm
    valid_min:  0
    valid_max:  20000
rhos_860 = ds["rhos"].sel({"wavelength_3d": 860}, method="nearest")

fig, ax = plt.subplots(figsize=(9, 5), subplot_kw={"projection": cartopy.crs.PlateCarree()})
ax.gridlines(draw_labels={"left": "y", "bottom": "x"}, linewidth=0.25)
ax.coastlines(linewidth=0.5)
ax.add_feature(cartopy.feature.OCEAN, edgecolor="w", linewidth=0.01)
ax.add_feature(cartopy.feature.LAND, edgecolor="w", linewidth=0.01)
rhos_860.plot(x="longitude", y="latitude", cmap="Greys_r", vmin=0, vmax=1.0)
plt.show()

rhos_860
<xarray.DataArray 'rhos' (number_of_lines: 1709, pixels_per_line: 1272)> Size: 9MB
[2173848 values with dtype=float32]
Coordinates:
    wavelength_3d  float64 8B 860.0
    longitude      (number_of_lines, pixels_per_line) float32 9MB ...
    latitude       (number_of_lines, pixels_per_line) float32 9MB ...
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:  Surface reflectance
    valid_min:  -0.05
    valid_max:  1.0
#Now flags 
ds["l2_flags"]
<xarray.DataArray 'l2_flags' (number_of_lines: 1709, pixels_per_line: 1272)> Size: 9MB
[2173848 values with dtype=int32]
Coordinates:
    longitude  (number_of_lines, pixels_per_line) float32 9MB ...
    latitude   (number_of_lines, pixels_per_line) float32 9MB ...
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:      Level-2 Processing Flags
    valid_min:      -2147483648
    valid_max:      2147483647
    flag_masks:     [          1           2           4           8         ...
    flag_meanings:  ATMFAIL LAND PRODWARN HIGLINT HILT HISATZEN COASTZ SPARE ...

Use the cf_xarray library to create a mask including only land and excluding clouds and ice.

# Create Mask
ds["l2_flags"].cf.is_flag_variable
cldwater_mask = (
    (ds["l2_flags"].cf == "LAND")
& ~(ds["l2_flags"].cf == "CLDICE")
)
# Apply mask, creating new dataset
rhos = ds["rhos"].where(cldwater_mask)
# Plot a band of masked data
rhos_860 = rhos.sel({"wavelength_3d": 860}, method="nearest")

fig, ax = plt.subplots(figsize=(9, 5), subplot_kw={"projection": cartopy.crs.PlateCarree()})
ax.coastlines(linewidth=0.25)
ax.gridlines(draw_labels={"left": "y", "bottom": "x"}, linewidth=0.25)
ax.add_feature(cartopy.feature.OCEAN, edgecolor="w", linewidth=0.01)
ax.add_feature(cartopy.feature.LAND, edgecolor="w", linewidth=0.01)
ax.add_feature(cartopy.feature.LAKES, edgecolor="k", linewidth=0.1)
rhos_860.plot(x="longitude", y="latitude", cmap="Greys_r", vmin=0, vmax=1.0)
plt.show()

To more easily work with the data we can project and resample it onto a regular grid.

To do this with rioxarray, we need to change the order of the dimensions in the array. Hyperspectral data is often stored (y,x,band) since most processing; however, rioxarray expects the dimensions to be in the order of (band, y, x), so we must transpose the array to grid the data using rioxarray.

# Transpose data
sr_src = rhos.transpose("wavelength_3d", ...)
sr_src
<xarray.DataArray 'rhos' (wavelength_3d: 122, number_of_lines: 1709,
                          pixels_per_line: 1272)> Size: 1GB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * wavelength_3d  (wavelength_3d) float64 976B 346.0 351.0 ... 2.258e+03
    longitude      (number_of_lines, pixels_per_line) float32 9MB -99.13 ... ...
    latitude       (number_of_lines, pixels_per_line) float32 9MB 13.6 ... 36.6
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:  Surface reflectance
    valid_min:  -0.05
    valid_max:  1.0
# # Add coordinates for longitude and latitude
# sr_src.coords["longitude"] = ds["longitude"]
# sr_src.coords["latitude"] = ds["latitude"]
       
sr_src = sr_src.rio.set_spatial_dims("pixels_per_line", "number_of_lines").rio.write_crs("epsg:4326")
sr_src
<xarray.DataArray 'rhos' (wavelength_3d: 122, number_of_lines: 1709,
                          pixels_per_line: 1272)> Size: 1GB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * wavelength_3d  (wavelength_3d) float64 976B 346.0 351.0 ... 2.258e+03
    longitude      (number_of_lines, pixels_per_line) float32 9MB -99.13 ... ...
    latitude       (number_of_lines, pixels_per_line) float32 9MB 13.6 ... 36.6
    spatial_ref    int64 8B 0
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:  Surface reflectance
    valid_min:  -0.05
    valid_max:  1.0
#If too many pixels too much memory could be used...........................

sr_dst = sr_src.rio.reproject(
    dst_crs=sr_src.rio.crs,
    src_geoloc_array=(
        sr_src.coords["longitude"],
        sr_src.coords["latitude"],
    ),
    nodata=np.nan
).rename({'y': 'latitude', 'x': 'longitude'})
sr_dst
<xarray.DataArray 'rhos' (wavelength_3d: 122, latitude: 1524, longitude: 2037)> Size: 2GB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * longitude      (longitude) float64 16kB -105.2 -105.2 ... -74.46 -74.44
  * latitude       (latitude) float64 12kB 36.62 36.6 36.59 ... 13.63 13.62 13.6
  * wavelength_3d  (wavelength_3d) float64 976B 346.0 351.0 ... 2.258e+03
    spatial_ref    int64 8B 0
Attributes:
    long_name:   Surface reflectance
    valid_min:   -0.05
    valid_max:   1.0
    _FillValue:  nan
def create_pace_rgb(rhos_ds, red=610, green=555, blue=465, pct_low=0.02, pct_high=0.98, gamma=0.5):

    
    ds_rgb = sr_dst.sel(wavelength_3d=[red,green, blue], method='nearest')

    # Create NaN Mask
    nan_mask = np.isnan(ds_rgb.data)

    # Stretch and Clip
    for i in range(ds_rgb.data.shape[0]): 
        channel = ds_rgb.data[i]
        valid = channel[~nan_mask[i]]
        if len(valid) > 0:
            p_low = np.percentile(valid, pct_low)
            p_high = np.percentile(valid, pct_high)
            channel = (channel - p_low) / (p_high - p_low)
            np.clip(channel, 0, 1, out=channel)
            
    # Apply Gamma and Clip
    ds_rgb.data = np.power(ds_rgb.data, 0.5)
    np.clip(ds_rgb.data, 0, 1, out=ds_rgb.data)
    return ds_rgb
# # Plot RGB Image - Still needs some brightening
# red_approx = 680.0
# green_approx = 540.0
# blue_approx = 450.0
rgb_da = create_pace_rgb(sr_dst, pct_low=0.01, pct_high=0.99, gamma=0.3)
rgb = rgb_da.hvplot.rgb(x='longitude', y='latitude', bands='wavelength_3d', aspect = 'equal', frame_height=400, geo=True, crs='EPSG:4326', title="RGB PACE Image with Pre-selected bands")
rgb
C:\Users\ebolch\AppData\Local\Temp\1\ipykernel_14144\2150097913.py:20: RuntimeWarning: invalid value encountered in power
  ds_rgb.data = np.power(ds_rgb.data, 0.5)
import holoviews as hv
import hvplot.xarray  # noqa
import geoviews as gv
import panel as pn
import pandas as pd
import xarray as xr
import hvplot.pandas

hv.extension('bokeh')

# Set Point Limit
POINT_LIMIT = 10

# Set up Color Cycling
color_cycle = hv.Cycle('Category20')

# Hack the colors to match - not sure why index adjustment is necessary
colors = color_cycle.values[1:POINT_LIMIT+1]

# Starting from known point with data (in EPSG:4326)
lon_mid = -85.3879
lat_mid = 31.3446

# Initialize points in EPSG:4326
initial_data = {'longitude': [lon_mid], 'latitude': [lat_mid]}
points = hv.Points(initial_data, kdims=['longitude', 'latitude'])

points_stream = hv.streams.PointDraw(
    data=points.columns(),
    source=points,
    drag=True,
    num_objects=POINT_LIMIT,
    styles={'fill_color': colors, 'line_color': 'gray'}
)

# Add vector coastlines (and optionally borders/land)
coastlines = gv.feature.coastline(projection=cartopy.crs.PlateCarree()).opts(line_color='black', line_width=1)

# Streams for hover and tap (source is the RGB image)
posxy = hv.streams.PointerXY(source=rgb, x=lon_mid, y=lat_mid)
tapxy = hv.streams.Tap(source=rgb, x=lon_mid, y=lat_mid)

# Subscriber to add point on tap (no transformation needed)
def add_point_on_tap(x, y):
    if x is None or y is None:
        return
    current_data = points_stream.data.copy()
    num_points = len(current_data['longitude'])
    if num_points >= POINT_LIMIT:
        return
    current_data['longitude'].append(x)
    current_data['latitude'].append(y)
    points_stream.send(data=current_data)

tapxy.add_subscriber(add_point_on_tap)

# Function to build spectral plots (vectorized; no transformation needed)
def spectra_cb(data, x, y):
    plots = []
    if data is None or not any(len(d) for d in data.values()):
        return hv.Curve([]).opts(
            xlabel='wavelength (nm)', ylabel='rhos', 
            frame_width=400, frame_height=500, title='Spectra'
        )
    
    # Points are in 4326
    lons = xr.DataArray(data['longitude'], dims='point')
    lats = xr.DataArray(data['latitude'], dims='point')
    selected = sr_dst.sel(longitude=lons, latitude=lats, method='nearest')
    selected = selected.assign_coords(point=[f'Point {i+1}' for i in range(selected.point.size)])
    
    plots.append(
        selected.hvplot.line(
            y='rhos', x='wavelength_3d', by='point', color=color_cycle
        )
    )
    
    # Add hover spectrum if position available
    if x is not None and y is not None:
        hover_data = sr_dst.sel(longitude=x, latitude=y, method="nearest")
        plots.append(
            hover_data.hvplot.line(
                y='rhos', x='wavelength_3d', 
                color='black', line_dash='dashed', label="Hover"
            )
        )
    
    return hv.Overlay(plots).opts(
        legend_position='top_left', show_legend=True
    )

# Define the Dynamic Map (updates on point_draw, hover)
spectra_dmap = hv.DynamicMap(spectra_cb, streams=[posxy, points_stream])

# Dynamic table for displaying points (in 4326)
def points_table(data):
    if data is None or not len(data['longitude']):
        return hv.Table(pd.DataFrame(columns=['ID', 'longitude', 'latitude'])).opts(width=300, height=200)
    df = pd.DataFrame({
        'ID': [f'Point {i+1}' for i in range(len(data['longitude']))],
        'longitude': data['longitude'], 
        'latitude': data['latitude']
    })
    return df.hvplot.table(title='Collected Points', width=300, height=200)

points_table_dmap = hv.DynamicMap(points_table, streams=[points_stream])

# Apply opts to points
points.opts(
    size=10, tools=['hover'], active_tools=['point_draw'], 
    color='white', line_color='gray'
)

# Button to download points as CSV (points in 4326)
def save_points():
    if not points_stream.data['longitude']:
        return None
    df = pd.DataFrame({
        'longitude': points_stream.data['longitude'], 
        'latitude': points_stream.data['latitude']
    })
    return df.to_csv(index=False).encode('utf-8')

download_button = pn.widgets.FileDownload(
    callback=save_points, 
    filename='collected_points.csv', 
    label='Download Points as CSV',
    button_type='primary'
)

# Layout: RGB + coastlines on left, spectra on right, table and button below
layout = pn.Column(
    pn.Row(rgb * coastlines * points, spectra_dmap),
    pn.Row(points_table_dmap, download_button)
)

# If using Jupyter, use layout.servable() or just return layout
layout
hv.extension('bokeh')

POINT_LIMIT = 10
color_cycle = hv.Cycle('Category20')

# Set up a holoviews points array to enable plotting of the clicked points
xmid = sr_dst.longitude.values[int(len(sr_dst.longitude) / 2)]
ymid = sr_dst.latitude.values[int(len(sr_dst.latitude) / 2)]

first_point = ([xmid], [ymid], [0])
points = hv.Points(first_point, vdims='id')
points_stream = hv.streams.PointDraw(
    data=points.columns(),
    source=points,
    drag=True,
    num_objects=POINT_LIMIT,
    styles={'fill_color': color_cycle.values[1:POINT_LIMIT+1], 'line_color': 'gray'}
)

posxy = hv.streams.PointerXY(source=rgb_map, x=xmid, y=ymid)
clickxy = hv.streams.Tap(source=rgb_map, x=xmid, y=ymid)

# Function to build spectral plot of clicked location to show on hover stream plot
def click_spectra(data):
    coordinates = []
    if data is None or not any(len(d) for d in data.values()):
        coordinates.append(clicked_points[0][0], clicked_points[1][0])
    else:
        coordinates = [c for c in zip(data['x'], data['y'])]
    
    plots = []
    for i, coords in enumerate(coordinates):
        x, y = coords
        data = sr_dst.sel(longitude=x, latitude=y, method="nearest")
        plots.append(
            data.hvplot.line(
                y="rhos",
                x="wavelength_3d",
                color=color_cycle,
                label=f"{i}"
            )
        )
        points_stream.data["id"][i] = i
    return hv.Overlay(plots)

def hover_spectra(x,y):
    return sr_dst.sel(longitude=x,latitude=y,method='nearest').hvplot.line(y='rhos',x='wavelength_3d',
                                                                            color='black', frame_width=400)
# Define the Dynamic Maps
click_dmap = hv.DynamicMap(click_spectra, streams=[points_stream])
hover_dmap = hv.DynamicMap(hover_spectra, streams=[posxy])
# Plot the Map and Dynamic Map side by side
hv.Layout(hover_dmap*click_dmap + rgb_map * points).cols(2).opts(
    hv.opts.Points(active_tools=['point_draw'], size=10, tools=['hover'], color='white', line_color='gray'),
    hv.opts.Overlay(show_legend=False, show_title=False, fontscale=1.5, frame_height=480)
)
# Plotting only to collect data points in points_stream
hv.extension('bokeh')

# Set Point Limit
POINT_LIMIT = 10

# Set up Color Cycling
color_cycle = hv.Cycle('Category20')
colors = [color_cycle[i] for i in range(5)]

#starting from known point with data 
xmid = -85.3879
ymid = 31.3446

# Initialize points (xmid, ymid are correct for sr_dst, in EPSG:4326)
first_point = ([xmid], [ymid], [0])
points = hv.Points(first_point, kdims=['x', 'y'], vdims='id')
points = points#.opts(projection=cartopy.crs.PlateCarree())  #setting this to this projection makes the background map dissapear. 


points_stream = hv.streams.PointDraw(
    data=points.columns(),
    source=points,
    drag=True,
    num_objects=POINT_LIMIT,
    styles={'fill_color': color_cycle.values[1:POINT_LIMIT+1], 'line_color': 'gray'},
    add=True
)

posxy = hv.streams.PointerXY(source=map, x=xmid, y=ymid) #change source
clickxy = hv.streams.Tap(source=map, x=xmid, y=ymid) #change source

# Function to build spectral plot of clicked location to show on hover stream plot
def click_spectra(data):
    coordinates = []
    if data is None or not any(len(d) for d in data.values()):
        coordinates.append((xmid, ymid))
    else:
        coordinates = [c for c in zip(data['x'], data['y'])]
        ids = [c for c in zip(data['id'])]
    
    plots = []
    for i, coords in enumerate(coordinates):
        x, y = coords
        idval = i
        data = sr_dst.sel(x=x, y=y, method="nearest")
        plots.append(
            data.hvplot.line(
                y='rhos',
                x='wavelength_3d',
                color=color_cycle,
                label=f"{i}"
            )
        )
        
        # if len(points_stream.data["id"] != num_points:
        #     print(f"Warning: Length mismatch. x: {num_points}, id: {len(id_data)}. Reassigning ids.")
        #     data['id'] = list(range(num_points))  # Assign incremental ids
       
        #if 'id' in idval < len(points_stream.data["id"]):
        points_stream.data['id'] = i
    return hv.Overlay()
    

# Define the Dynamic Maps
click_dmap = hv.DynamicMap(click_spectra, streams=[points_stream])


# Plot the Map and Dynamic Map side by side
layout = hv.Layout(click_dmap + map * points ).opts(
    hv.opts.Points(active_tools=['point_draw'], size=10, tools=['hover'], color='white', line_color='gray'),
    hv.opts.Overlay(show_legend=False, show_title=False, fontscale=1.5, frame_height=480)
)

# Debug saved coordinates
print(f"Initial points_stream.data: {points_stream.data}")

# Display the layout
layout
print(points_stream)
# # Define transformer for EPSG:3857 (map clicks) to EPSG:4326 (sr_dst)
transformer = pyproj.Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)

# # Function to check if coordinates are in degree range (EPSG:4326)
def is_degree_range(x, y):
    return -180 <= x <= 180 and -90 <= y <= 90

plots = []
for i, (x, y, id_val, color) in enumerate(zip(points_stream.data['x'], points_stream.data['y'], points_stream.data['id'], points_stream.data['fill_color'])):
    print(f"Point {i}: x={x}, y={y}, id={id_val}, color={color}")
    if is_degree_range(x, y): 
        data = sr_dst.sel(x=x, y=y, method="nearest")
        plots.append(data.hvplot.line(
                    y='rhos',
                    x='wavelength_3d',
                    color=color,
                    label=f"{id_val}"
                ))
    else:
                # Transform other points from EPSG:3857 to EPSG:4326
                x_trans, y_trans = transformer.transform(x, y)
                print(f"Transformed ({x}, {y}) to ({x_trans}, {y_trans})") 
                data = sr_dst.sel(x=x_trans, y=y_trans, method="nearest")
                plots.append(data.hvplot.line(
                    y='rhos',
                    x='wavelength_3d',
                    color=color,
                    label=f"{id_val}"
                ))
    
hv.Overlay(plots) if plots else hv.Overlay([])

Calculate unique indices

Hyperspectral-enabled Vegetation Indices

We’ll take the Chlorophyll Index Red Edge (CIRE) as an example of a hyperspectral-enabled VI. CIRE uses bands from the red edge and the NIR to get at relative canopy chlorophyll content.

CIRE: (ρ 800 /ρ 705)−1

Carotenoid Content Index (Car): [(1/ρ495)−(1/ρ705)] * ρ800

Information on Hyperspectral enabled indices by PACE here: https://www.earthdata.nasa.gov/apt/documents/landvi/v1.0#mathematical_theory

## Chlorophyll Index Red Edge (CIRE)
# Select bands
sr_800 = sr_dst.sel({"wavelength_3d": 800}, method="nearest")
sr_705 = sr_dst.sel({"wavelength_3d": 705}, method="nearest")
#Calculate
cire = (sr_800 / sr_705) - 1
cire.attrs["long_name"] = "Chlorophyll Index Red Edge (CIRE)"
#Plot
map1 = cire.hvplot.image(x='x', y='y', frame_height=400, geo=True, cmap="viridis", tiles='CartoDark', title="CIRE from PACE OCI")
map1
## Carotenoid Content Index (Car) 
# Select additional bands
sr_495 = sr_dst.sel({"wavelength_3d": 495}, method="nearest")

#Calculate
car = ((1 / sr_495)- (1 / sr_705)) * sr_800
car.attrs["long_name"] = "Caretenoid Content Index (Car)"
map2 = car.hvplot.image(x='x', y='y', frame_height=400, geo=True, cmap="plasma", tiles='CartoDark', title="CAR from PACE OCI")
map2
# TODO No need to use hv layout for side-by side or layered on top, can simplify to this and providing map_opts will suppress warning.
map1.opts(frame_height=400, frame_width=480)+map2.opts(frame_height=400, frame_width=480)
# layout = hv.Layout(map1 + map2 ).cols(2).opts(hv.opts.Image(width=500, height=480))
# layout
#Enjoy PACE data 
import holoviews.element.tiles as hv_tiles
print(list(hv_tiles.tile_sources.keys()))